Levy distribution in many-particle quantum systems 
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Levy distribution, previously used to describe complex behavior of classical systems, is shown to 
characterize that of quantum many-body systems. Using two complimentary approaches, the canon- 
ical and grand-canonical formalisms, we discovered that the momentum profile of a Tonks-Girardeau 
gas, - a one-dimensional gas of N impenetrable (hard-core) bosons, harmonically confined on a lat- 
tice at finite temperatures, obeys Levy distribution. Finally, we extend our analysis to different 
confinement setups and demonstrate that the tunable Levy distribution properly reproduces mo- 
mentum profiles in experimentally accessible regions. Our finding allows for calibration of complex 
many-body quantum states by using a unique scaling exponent. 

PACS numbers: 05.40. Fb; 67.85.-d; 47.27.eb; 05.30.Jp 



I. INTRODUCTION 

Since the first observations of quantum collective phe- 
nomena, the quantum systems with strongly interacting 
constituents have become of paramount interest in con- 
dense matter community A new wave of activity 
in the area of quantum many-body systems has been 
burgeoning with the advent of laser cooling techniques 
0. Many of those quantum models, which were thought 
to be theoretical abstractions, have since been imple- 
mented with cold atoms [!, One of such models, a 
one-dimensional system of hard-core bosons - a Tonks- 
Girardeau (TG) gas Ja|_- has been thoroughly probed in 
recent experiments [fj, \1\ . 

While the density profile and the energy spectrum 
of hard-core bosons resemble that of non-interacting 
fcrmions, its momentum distribution (MD) exhibits dis- 
tinct features. The ground state of a homogeneous TG 
gas is known to possess an infrared divergence in the 
thermodynamic limit, N — > oo, n(p) oc p~ 1 / 2 [8j, which, 
however, vanishes upon addition of a harmonic confine- 
ment 0. So far, there are no analytic results on the in- 
frared behavior of the finitc-tcmperature MD for a finite 
number of bosons in a harmonic trap. In a sole harmonic 
confinement, the ground-state MD decays as n(p) oc p~ 4 
at the high momentum regime fioj . Yet in the presence of 
an optical lattice, which sets an upper momentum scale 
given by the recoil momentum, hk^, where fcz, is the wave 
vector of the laser beam, this region cannot be resolved 
in present state-of-the-art experiments. 

The finite-temperature MDs of TV < 20 strongly repul- 
sive bosons, confined on a one-dimensional optical lat- 
tice and an additional harmonic trap (a so called "Id 
tube"), have been measured experimentally [6|. These 
measurements revealed that in the intermediate region, 
0-3 < \p/Hki,\ < 1, a momentum profile can be approx- 
imated by a power-law, n(p) oc p~ 7 . The exponent 7 
depends on temperature, density, and strength of the 
atom- atom interactions. Results of numerical Monte- 
Carlo simulations have corroborated the experimental 
finding [Til ]. 

What kind of momentum distribution emerges in a sys- 



tem of strongly repulsive bosons? In the present work, 
we attest that Levy distribution describe the MD of N 
thermalized hard-core bosons in various one-dimensional 
confinements, in particular, within a single Id tube as 
well as within array of Id tubes probed experimentally 
[1, 0] • Levy statistics [l2| are known to describe classical 
chaotic transport [H[ , processes of subrecoil laser cooling 
[l4| , fluctuations of stock market indices [l5[ , time series 
of single molecule blinking events |16[ , or bursting activ- 
ity of small neuronal networks [17] . The appearance of 
Levy distribution in a system output is a strong indica- 
tor of a long-range correlation "skeleton" which conducts 
system intrinsic dynamics fPH , Il8j . However the Levy dis- 
tribution has at no time emerged in the context of many- 
particle quantum systems before. The great advantage of 
the Levy-based analysis is its capability of calibration of 
the TG in different quantum regimes by a unique scaling 
exponent a. 



The paper is organized as follows. In Sec. II, we intro- 
duce the lattice TG model and employ two complimen- 
tary approaches to study its finite-temperature proper- 
ties. In Sec. Ill, we fit exact results of the calculation 
within the grand-canonical formalism [l9j by the Levy 
distribution. To provide an insight into experimental sit- 
uation we devoted Sec. IV to the analysis of momen- 
tum profiles averaged over the array of Id tubes. We 
proceed with the results obtained within the canonical 
formalism [20] | and conclude the section with the analy- 
sis of the experimental data from the Ref. [6j. In Sec. V, 
we demonstrate the universality of the Levy spline ap- 
proximation by addressing the MD of a TG gas in one- 
dimensional confinements of various geometry. We elab- 
orate on the case of a TG gas in a sole harmonic con- 
finement, in a box, and in a sole optical lattice with im- 
penetrable (hard wall) boundaries. Finally, in Sec. VI, 
we summarize our results. Some of important technical 
details are deferred to the Appendixes A and B. 
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FIG. 1: (Color online) Normalized momentum density profiles 
(MD) for N = 40(a); 30(b); 10(c); 90(d) hard-core bosons in 
a Id tube, thin (black) lines, at four different sets of param- 
eters are compared to Levy distributions, thick (red) lines. 
Dashed-doted line (c) corresponds to a Gaussian approxima- 
tion. The parameters (N/n s ,kBT/ J) are indicated by sym- 
bols T, ■. • in Fig. 2. The thin solid (blue) icon indicates 
the confinement geometry. 



To find the reduced single-particle density matrix, we 
employ here the grand-canonical [19| and the canonical 
formalism [5o| (See Appendixes A, B). The first is rele- 
vant for a system being in contact with a thermal cloud 
at constant temperature T, while the second describes an 
isolated many-particle system. 

The difference between the momentum distributions 
obtained within the grand-canonical and the canonical 
descriptions is mediated by the number of particles and 
becomes negligible for TV > 10 [l9j]. In addition, at finite 
temperatures fcgT > 0.1J, systems with different num- 
ber of particles, but the same densities N/n s , possess the 
same momentum profiles Here, n s — 8(J^) 1 / 2 /7r is 
the number of single-particle eigenstates with non-zero 
population at the trap center (22|. The latter also yields 
the critical number of bosons in a Id tube required to 
form the Mott-insulator in the trap center at zero tem- 
perature. 

The typical examples of the MD obtained here within 
the grand-canonical formalism are presented in Fig. 1 for 
N > 10, where the results of both the grand-canonical 
and canonical descriptions practically identical. 



II. TONKS-GIRARDEAU GAS AT FINITE 
TEMPERATURES ON A ID OPTICAL LATTICE 

A Bose gas confined in a deep optical lattice is well 
described by the Bose-Hubbard Hamiltonian [2l| 

H = -jY^{b\h+i + h.c) + uY^l 2 ni + -^n l {n l -l), 
i i i 

(1) 

where ni = bjbi is the particle number operator on lat- 
tice site I, J is the hopping strength, the parabolicity 
v = AIuqcP 1 2 is the amplitude of the external harmonic 
potential with the trapping frequency uiq, M is the mass 
of the atom, and d is the lattice constant. The last term 
in ([1]) describes the on-site atom-atom interactions. Here 
we are interested in the TG regime, where the strength 
of the repulsive atom-atom interaction considerably ex- 
ceeds the kinetic (hopping) energy, i.e., U/J — > oo, @. 
Therefore, the interaction term can be substituted by the 
condition that two bosons cannot occupy the same lattice 
site 

The momentum distribution can be obtained from the 
reduced single-particle density matrix p n ,i, reading 

^) = i*(p)i 3 E e "* (B "Vi. ( 2 ) 

where the momentum p is written in units of recoil mo- 
mentum hkL = nh/d, and the envelope 3>(p) is the 
Fourier transform of the Wannier function. Note that 
the latter, as well as the hopping J, is solely denned by 
the lattice depth, Vb, measured in units of recoil energy, 
Er = (hk L ) 2 /2M. 



III. LEVY SPLINE APPROXIMATION 

In view of the strong non-Gaussian behavior of the 
function n(p) (see Fig. lc) and an apparent power-law in- 
termediate asymptotics [y, [ll| , it is tempting to compare 
the MD of a TG gas with symmetric Levy distribution 
[T^ |. The latter, L a (p), is a natural generalization of the 
Gaussian distribution, and it is defined by the Fourier 
transform of its characteristic function; i.e., 

1 f°° 

L a ( P ) = —J HO^dt, = exp(-C a \tn , 

where the exponent < a < 2, and C a > is some con- 
stant [12] . For a — 2 the distribution is the Gaussian. 
The value of a is invariant under the scaling transforma- 
tion p — > const ■ p, which can be used for fitting. 

The Levy distribution exhibits a power-law asymp- 
totics, L a (p) oc p~( a+1 \ as \p\ — > oo. These "heavy" 
tails cause the variance of Levy distributions to diverge 
for all a < 2. However, this asymptotic limit is not rel- 
evant for our objective: the MD for the system (1) is 
bounded by the width of |$(p)|, (V / 'Er^Ncl [ill, so 
that the MD variance remains finite. 

The fitting of the calculated MD n{p) by the Levy 
distribution L a (p) was automatized and performed nu- 
merically by the global minimum search of the root mean 

square deviation, A = {X)„[n(p n ) - L a (p n )} 2 } 1/2 /n(0), 
in a 2d parametric space {a, C a } ([3]). For the results 
shown on Figs. 1,2, we took 50 spline points equally 
spaced on the interval p € (Ojhk^). The initial area, 
a e (0, 2), C a € (0, 30), covered by 10 x 10 grid was iter- 
atively converged to the global minimum by decreasing 
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FIG. 2: (Color online) The scaling exponent of Levy distribu- 
tion, a (top), and the corresponding root mean square devi- 
ation from the exact momentum density profile, A (bottom) , 
of TG gas as a function of temperature and particle density 
for two different amplitudes of the optical lattice, Vo = 4.6Er 
(left column), and Vo = 9.3Er (right column), see text for 
more details. The amplitude of harmonic confinement is cho- 
sen as v — 2 ■ I0 -4 -Er. 



the grid step in a and C a until the desired relative accu- 
racy (fixed to 1% in all figures) was reached. 

In Fig. 2, we show the dependence of scaling exponent 
a (top) on the scaled temperature ksT/J and the scaled 
particle density, N/n s (obtained for N = 10, 100, and 
n s = 48.78), together with the associated mean square 
deviation A (bottom), at relatively low and high am- 
plitude of the optical lattice. These diagrams (see also 
Fig. 1, for concrete examples of the MD and their Levy- 
spline approximations) constitute the first main result of 
this work, namely, the convergence of the MD of a TG 
gas towards the Levy distribution with increasing the 
temperature. 

Despite of infinitely strong on-site repulsive interac- 
tions, the systematic increase of the Levy exponent a 
with temperature, cf. Fig. 2 (top) , is consistent with the 
high-temperature limit where the ideal Bose gas obeys 
the classical Boltzmann-Maxwell statistics with a Gaus- 
sian MD, i.e., a = 2. 



IV. 



AVERAGING OVER THE ARRAY OF ID 
OPTICAL LATTICES 



The single Id tube realization, using a fixed temper- 
ature T and a fixed particle number N, is not directly 
accessible with the present state-of-art experiments. The 
averaging over an array of tubes done in Ref. @ can be 
understood as an averaging over many realizations, with 
different parameters, T and N. 

Namely, the experimental setup [6] produces an ar- 
ray of independent Id tubes, with different numbers of 
particles, iV,-. The probability of having a tube with TV 
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FIG. 3: (Color online) Unnormalized momentum distribu- 
tions (MD) for N = 3,..., 18 hard-core bosons, thin (black) 
lines, obtained within the canonical approach (|B1I) and the 
MD averaged over the array of Id tubes according to the ex- 
perimental procedure used in Ref. p\, thick (red) line, with 
its best Levy-spline, thick (blue) dashed line. 



particles is given by @ 
g(N) 1 



3/V c 2/3 iVV3' 



N < N r 



(4) 



where the number of particles in a central tube N c is a 
unique parameter. Assuming the same initial temper- 
ature in all the tubes for a shallow Id lattice potential, 
Vi n , the temperatures at the experimentally adjusted lat- 
tice depth, Vq , can be obtained by using the conservation 
of entropy in each tube during the subsequent adiabatic 
increase of the lattice depth from V ln to Vq. Therefore, 
tubes with the different number of particles, N%, acquire 
different final temperatures, TJ, at V > V[ n ((J. 

We implemented this averaging procedure with N c — 
18, for a set of individual MDs pre-calculated within the 
canonical formalism. The result for fc^T = 0.5J, Vq — 
V in = 4.6£ fl and v = 8 • lO" 4 ^ is depicted in Fig. 3. 

Surprisingly enough, the averaged MD can be perfectly 
approximated by the Levy distribution with an "average" 
scaling exponent, a; this is so despite the sizable disper- 
sion of a-values appearing in the different momentum 
profiles. 

Therefore, the averaged momentum profile appears as 
a superposition of several different profiles. It is known 
that using a proper weight function, /(/?), one can con- 
struct a Levy distribution from a parameterized set of 
Gaussian distributions of different dispersion, L a (p) = 
X!^ G(p,p)f(fi)dp [13. Yet, to the best of our knowl- 
edge, there are no results concerning the superposition of 
many different Levy distribution functions with different 
exponents a. 

A next objective is the comparison of our scheme with 
the experimental data from Ref. Q, see Fig. 4. As one 
can deduce the Levy distribution yields an excellent ap- 
proximation for the MD of the experimental system al- 
though the latter does not map precisely to a TG gas, 
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FIG. 4: (Color online) The experimental data for the MD 
of the Id quantum gases for different axial lattice depths [f| , 
grey dots; Levy distributions with the scaling exponent a, 
solid (red) lines. The thin solid (blue) icon indicates the con- 
finement geometry. 



but rather corresponds to a set of N soft-core bosons 
with U/J < oo (HJ) 11]. In the experiment, one used a 



on an intermediate range 



1 = 1.2; 0.9; 0.4; 0.2 



power-law fit, n(j>) oc p 
(0.3 < \ P /hk L \ < 1) [6], with 7cxp - 
for the data shown on Figs. 4(a-d), respectively. We em- 
phasize, however, that in the experimentally accessible 
region, \p\ < hk^, the power law behavior with the ex- 
ponent 7 CX p — 1 = a of the Levy distribution is not yet 
valid, but is assumed only for much larger momentum 
values. Therefore, the theoretical estimates a in Fig. 4 
exceed those intermediate range power-law fit- values, i.e. 
one consistently finds that a > 7 CX p — 1- 



V. UNIVERSALITY OF THE LEVY-SPLINE 
FOR A TONKS-GIRARDEAU GAS IN ID 
CONFINEMENT POTENTIALS 

So far we have been dealing with a TG gas on a lattice 
with an additional harmonic potential (fT]), which is the 
only confinement where the MD of a TG gas was experi- 
mentally studied in detail [6]. In the present section, we 
demonstrate applicability of the Levy-spline approxima- 
tion to the numerically obtained MD of hard-core bosons 
at finite temperatures in various Id confinements: a sole 
harmonic trap, Fig. [SJ a box, Fig. (Ha,b), and a sole op- 
tical lattice with impenetrable boundaries, Fig. [BJc,d). 

A TG gas confined in a general potential, V(x), is de- 
scribed by the sum of the single particle Hamiltonians, 



H 



N 

£ 

i=l 



n 2 d 2 

'2M dx 2 



(5) 



with the hard-core constrain on the bosonic many- 
particle wave function: ^(xx, X2, xn) = if \xi — Xj < 
a, where Xi is the position of ith particle, and a is the 
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FIG. 5: (Color online) Normalized momentum density pro- 
files (MD) for TV = 5(a,b); 10(c,d) hard-core bosons in a 
sole id harmonic trap, thin (black) lines, at the two different 
temperatures, ksT '/Hljq = 1 (a,c), ksT/trwo — 1.5 (b,d) are 
compared to Levy distributions, thick (red) lines. Dashed line 
(c) underlines well known universal power-law asymptotic be- 
havior of the Tonks- Girardeau gas momentum distribution at 
large momenta, p~ 4 , while dashed-dotted line illustrates fail- 
ure of a Gaussian approximation. The thin solid (blue) icon 
indicates the confinement geometry. 



Id hard-core diameter. At low density, hard-bosons can 
be approximated by impenetrable point size particles, so 
that a = 0|. 

The MD of N hard-bosons in a box, Vbox(^i) = if < 
Xi < Lbox an d VboK — 00 otherwise, and in a harmonic 
confinement, V osc — mui 2 x 2 /2, was obtained numerically 
within the canonical formalism with the help of efficient 
method [25[ , which is detailed in Appendix IB 21 In the 
case of a sole lattice potential with L lattice sites, the MD 
was calculated within the grand-canonical formalism by 
using the same algorithm that was employed in Sec. IV 
(see Appendix A). 

It is known, that the MD of a TG gas in a homogeneous 
toroidal trap and in a harmonic confinement (Fig. [SJ ex- 
hibits a power-law behavior, p~ 4 , for p — > oo 8]. We 
also found the same power-law behavior for a TG gas in 
a box, see Fig. |UJa,b). This in fact turns out to be a 
general feature of a TG, which also persists at finite tem- 
peratures. However, the power-law tail, p~ a ~ l , of the 
Levy distribution with exponent a cannot decay faster 
than p~ 3 [l2j]. Therefore, our Levy-spline approxima- 
tion is exclusively aimed to a finite momentum region 
(0,p c ), where the asymptotic p~ 4 behavior has not yet 
developed. With Figs. [SJB1 we demonstrate that the MD 
of a TG gas at finite temperatures in various Id con- 
finement potentials (thin lines) can be perfectly approx- 
imated by Levy-splines (thick lines) over significant mo- 
mentum range. 

While at zero temperature we detect a notable devia- 
tion of the Levy-spline from the actual MD, with increas- 
ing the temperature the deviation becomes practically 
negligible. The exponent of Levy-splines for a TG gas 
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FIG. 6: (Color online) Normalized momentum density pro- 
files (MD) for TV = 10(a,d); 5(b); 50(c) hard-core bosons 
in a Id box (a,b) of size I/box an d on a lattice (c,d) with 
L — 100 lattice sites and the potential amplitudes Vo = 
4.6Er(c); 9.9Eu(d), thin (black) lines, at the temperatures 
fcsT/^box = 20 (a); 10 (b), fesT/J = 0.5 (c,d) are compared 
to Levy distributions, thick (red) lines. Here the frequency 
i^box = ft7r 2 /2MZ/box sets the energy scale, hwo, for a par- 
ticle in a box. The momentum, p, is measured in units of 
po — ft/Lbox(a,b) and po = HUl (c,d), in the case of a lattice 
potential and a box, respectively. The thin solid (blue) icons 
indicate the two corresponding types of confinement geome- 
try. 



We want to emphasize one aspect which we consider 
to be crucial for the understanding of our approach. It 
is known that the Levy distribution has appeared as the 
stable distribution, i.e. the "attractor" for normalized 
sums of independent and identically-distributed random 
variables, with no finite mean values. The latter has 
opened the door to anomalous statistics [T3 |. In our ap- 
proach, we employ the Levy distribution ([3]) merely as a 
mathematical function, which is a natural generalization 
of the standard Gaussian function, thus leaving aside all 
the possible statistical interpretations and speculations 
concerning physical processes responsible for the MD for- 
mation. A good example for such "applied" approaches is 
the widespread use of the Gaussian distribution, which 
describes the wave function of the ground state of the 
quantum harmonic oscillator [261 ] . the Green's function 
solution of the deterministic heat or diffusion equation 
[27| . and frequently results as a form function in many 
different physical contexts. 

While it is intuitively clear that it is the long-range 
correlations in the system that cause the emergence of 
Levy distributions, the task to unravel the inherent phys- 
ical mechanism(s) yielding this anomalous distribution 
remains a challenge. We think that the analysis of a 
reduced single-particle density matrix in a spirit of the 
theory of random Levy matrices [28[ may shed additional 
light on this intriguing issue. 



in the confinement potentials presented here depends, in 
general, on the confinement characteristics, on the gas 
temperature and on the particle number or the filling 
factor, N/L (in the case of a lattice). A more detailed 
analysis of this dependence will be addressed elsewhere. 

VI. CONCLUSIONS 

We have presented a study of the finite-temperature 
momentum distribution (MD) of Tonks-Girardeau (TG) 
gases confined in various Id confinement potentials. The 
MD of a TG gas on an optical lattice with a superimposed 
harmonic confinement depends on the particle density, 
on the lattice depth and on the gas temperature. We 
have shown that the tunable Levy distribution fits mo- 
mentum profiles up to one recoil momentum with high 
accuracy. This allows for calibration of TG states with 
a unique scaling exponent. Thus, our approach com- 
pletes the attempts to quantify the finite-temperature 
MDs by using a power-law fitting for the intermediate 
region 0.3 < \p/hk L \ < 1 @,[Tl|. 

We demonstrate that the MD of a TG gas confined 
in a Id box, in a sole Id harmonic potential and on a 
Id optical lattice with hard wall boundaries can be ap- 
proximated by the Levy distribution function on a 
finite momentum region (0,p c ), where for a TG gas the 
universal power-law, p -4 , at large momenta has not yet 
developed. We thus conjecture that Levy scaling of MD 
is generic feature of a thermalized TG gas. 
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Appendix A: Grand-canonical formalism 

Within the grand-canonical description, the number 
of bosons, TV, is a fluctuating quantity and the reduced 
single-particle density matrix is defined as the trace over 
the Fock space, 

p,u = ±Tr{(blb l )e-K H ->X i n i )} , ( A l) 

where the chemical potential /x is fixed to give the re- 
quired number of particles N — pij in the system, 
/3 = 1/ksT with ks being the Boltzmann factor, and 
Z is the grand-canonical partition function. The lat- 
ter coincides with that of non- interacting fermions, i.e., 
Z = J\j[l + e - ^^ - ^], where Cj stands for the single- 
particle energy spectrum. 

The trace over the the Fock space (|A1[) can be eval- 
uated exactly [19] by mapping the problem of hard- 
core bosons on that of spinless fermions via the Jordan- 
Wigner transformation [29j . We next use the expressions 
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for the elements of the reduced single-particle density 
matrix elaborated in [13, i.e., 



Pn,i = \ {dct [I + (I + A(n, 0)D(rc, /)] 
-det[I + D(n,0]} ,iovnj=l, 

D(n,Z) = 0(Z)Ue-( E -^/' £BT U t O(n); 



and 



Pn : n 



U I 



,(E-^I)/fc B T 



ut 



(A2) 
(A3) 

(A4) 



for the main diagonal elements. All operators entering 
above are the square matrices N x N defined as fol- 
lows: I denotes the identity matrix, 0(n)i,j = {— 5i d , i < 
ra— 1; Sij, i > Ti—1}, A(n, — 6i tn Sj,i, U is the orthog- 
onal matrix of eigenvectors satisfying the eigenproblem, 
H{U = UE, with single-particle version of Hamiltonian 



Hi = -jJ2(\ 1 + + h - c -) + v ^ l ' z \ l )( l i 



(A5) 



'E' 2 

i i 
and E is diagonal matrix of its eigenvalues, i.e., E<j = 
EiSij. Thus, to obtain the entire matrix p n j (pi >n = p* n ; ) 
one has to compute N(N — l)/2 determinants of N x A 
matrices. 



Appendix B: Canonical formalism 

In the canonical formalism, the particle reservoir is 
absent, meaning that one has to find the number, A/" e , 
self-consistently from those many particle states, \^ m ), 
together with their eigenenergies, E m , that contribute to 
the thermal superposition at a given temperature. The 
reduced single-particle density matrix is then obtained as 
a sum of thermally weighted density matrices, p m , eval- 
uated for each from Af e eigenstates separately; i.e., 



rn—1 



(Bl) 



The canonical partition function Z is expressed through 
the true many particle energy spectrum: Z = 

ipMe -0E m 

The energy spectrum of hard-core bosons is the same 
as that for spinless fermions: E m = Y^n=i e ot m , where 
a™ denotes the numbers of single-particle eigenlevels oc- 
cupied in the m-th many-particle state. 



1. Hard-core bosons on a lattice 

In the case of a lattice confinement, p m = p™ ; = 
( i Sf m \bl l bi\^f m ). The corresponding many-particle eigen- 
states can be represented as 



i* ro )=nE u ^^i°)' n=i '---' iv 



(B2) 



where U is complete orthogonal set of single-particle 
eigenvectors (see Appendix A). 

Each contribution, p™ , is related to the Green func- 
tion, G" 



(B3) 



Using the Jordan- Wigner transformation [29], the 
bosonic Green function can be rewritten as scaler prod- 
uct of two fermionic wave-functions, and subsequently as 
a determinant of matrix product [20|: 

G™ =dct[(P m (i))tp m (j)] , (B4) 

where the matrix P m {i)i,n has N + 1 columns 

C -Uj, a m, fori < i - l,7i= 1,..., N; 
P m (»)i,>i = | U/, am , for/ > i - l,7i = 1,.. . , N; 
{ S it i, forn = 7V + l. 

(B5) 

In comparison to the grand-canonical approach, the num- 
ber of operations needed to obtain the entire matrix pi d 
is a factor of Af e larger, and is growing with increasing 
temperature. 



2. Hard-core bosons in arbitrary confinement 
potential 

Above we detailed the canonical formalism for hard- 
core bosons on a lattice. Here, instead, we elaborate on 
a general case with arbitrary confinement, V(x). In this 
case, the momentum distribution, n(p), can be obtained 
from the reduced single-particle density matrix in the 
continuum: n(p) — (2t:)~ 1 J dxdye~ lk ( x ~ y ^ p(x,y), where 
p(x, y) is given by the thermal superposition (|B1[) with 
p m = p rn (x,y) defined below. 

To calculate the reduced single-particle density ma- 
trix, p m (x, y), in the continuum we make use of efficient 
method [25|, which represents it in terms of the single 
particle states: 



P m (x,y) 



N 

E 

ij'=l 



r a Ax)A m {x,y) itj i> aT {y), (B6) 



where N x N matrix A m (x,y) associated to the m-th 
many-particle state is 



A m (x,y) id = (-l) i+i dBbQ m (i,j) 



(B7) 



Q m (i, j) is a minor of matrix Q m obtained by crossing 
i-th row and j-th column, and matrix Q m itself is given 
as 



{x,y)i 



dx'ip* am (x')ip aT (x'), (B8) 



where ip n {x) are the eigenfunctions of single particle 
eigenproblem in a given trapping potential V(x). In 
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(IB8p . it is assumed that x > y, while for i < y: 
p m (x,y) = p m (y,x)*. The set of a™, as before in Ap- 
pendix IB 1| denotes the numbers of single-particle eigen- 
levels occupied in the m-th many-particle state. Addi- 
tionally, whenever det Q m ^ 0, (|B7|) can be represented 

as A" 1 = (Q^ l 1 ) T det Q m [Ifl, which is more efficient 
when implemented numerically. 

For a box and for a harmonic confinement, 
ipn(%) = \/2/Lhax sia(rmx/Lbax) and ip n (x) = 
\J\j2 n n\ (muJo/Trh) 1 ^ 4 exp ( y —mujQX 2 /2fi) H n (^muJo/hx), 
respectively, where H n (y) are the Hermite polynomials. 



For these confinements, the integral in (IB8|) . and the 
density matrices p m (x,y) can be obtained analytically 
for a small number of particles with the help of symbolic 
computational routines. However, expanding out m 
times N 2 determinants of (N — 1) X (N — 1) matrices for 
N > 5 quickly becomes unwieldy. Therefore, in Sec. V, 
the reduced density matrices p m (x,y) were calculated 
on a numerical grid. The grid step was sufficiently 
small to ensure a smooth representation of all single 
particle wave-functions, ip m (x), participating in excited 
many-particle states at a given temperature. 
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